library(tidyverse)
library(lme4)
library(lmerTest)
library(knitr)
library(MuMIn)
library(cowplot)
library(caret)
library(ROCR)algorithm = c("#006989", "#FEC601", "#F43C13", "#00A5CF", "#00A878")
instruction = wesanderson::wes_palette("Darjeeling1", 2, "continuous")
craving = wesanderson::wes_palette("Darjeeling1", 3, "continuous")
rating = c("#00A08A", "#F2AD00", "#F98400", "#FF0000")
dc_bw = readRDS("~/dc_bw.Rds")source("~/Documents/code/sanlab/DEV_scripts/behavioral/descriptives/ROC/load_data.R")Check if wrong buttons were used (i.e., not 5-8)
subs = data.all %>%
group_by(subjectID, run, rating) %>%
summarize(n = n()) %>%
spread(rating, n) %>%
mutate(messed = ifelse(is.na(`5`) & !is.na(`<NA>`), "yes", NA)) %>%
filter(messed == "yes") %>%
ungroup() %>%
select(subjectID) %>%
unique()
data.all %>%
group_by(subjectID, run, rating) %>%
summarize(n = n()) %>%
spread(rating, n) %>%
mutate(messed = ifelse(is.na(`5`) & !is.na(`<NA>`), "yes", NA)) %>%
filter(subjectID %in% subs$subjectID)Recoding
* DEV069: recode runs1
data.ex = data.all %>%
mutate(rating = ifelse(subjectID == "DEV069" & run == "run1", rating - 1, rating),
rating = ifelse(subjectID == "DEV069" & run == "run1" & is.na(rating), 8, rating),
rating = rating - 4) %>%
group_by(subjectID, wave) %>%
arrange(subjectID, run) %>%
mutate(trial = row_number())striping = read.csv("~/Documents/code/sanlab/DEV_scripts/fMRI/betaseries/ROC/striping_QC.csv")file_dir = "~/Documents/code/sanlab/DEV_scripts/fMRI/betaseries/ROC/dotProducts_ROC/"
file_pattern = "DEV[0-9]{3}_meanIntensity.txt"
file_list = list.files(file_dir, pattern = file_pattern)
intensities = data.frame()
for (file in file_list) {
temp = tryCatch(read.table(file.path(file_dir,file), fill = TRUE) %>%
rename("subjectID" = V1,
"meanIntensity" = V3) %>%
extract(V2, "beta", "beta_([0-9]{4}).nii") %>%
mutate(beta = as.integer(beta)), error = function(e) message(file))
intensities = rbind(intensities, temp)
rm(temp)
}file_dir = "~/Documents/code/sanlab/DEV_scripts/fMRI/betaseries/ROC/dotProducts_ROC/"
file_pattern = "DEV[0-9]{3}_dotProducts.txt"
file_list = list.files(file_dir, pattern = file_pattern)
dots = data.frame()
for (file in file_list) {
temp = tryCatch(read.table(file.path(file_dir,file), fill = TRUE) %>%
rename("subjectID" = V1,
"map" = V3,
"dotProduct" = V4) %>%
extract(V2, "beta", "beta_([0-9]{4}).nii") %>%
extract(map, "algorithm", "(.*)_.*.nii") %>%
mutate(beta = as.integer(beta)), error = function(e) message(file))
dots = rbind(dots, temp)
rm(temp)
}dots.merged = dots %>%
left_join(., intensities, by = c("subjectID", "beta")) %>%
group_by(subjectID, algorithm) %>%
mutate(rownum = row_number())
# plot original
dots.merged %>%
filter(algorithm == "logistic") %>%
ggplot(aes(1, meanIntensity)) +
geom_boxplot()# assess extreme values and exclude when calculating SDs
dots.merged %>%
filter(algorithm == "logistic") %>%
arrange(meanIntensity)dots.merged %>%
filter(algorithm == "logistic") %>%
arrange(-meanIntensity)# recode outliers as NA
dots.merged = dots.merged %>%
ungroup() %>%
mutate(meanIntensity = ifelse(meanIntensity > 1 | meanIntensity < -1, NA, meanIntensity),
median = median(meanIntensity, na.rm = TRUE),
sd3 = 3*sd(meanIntensity, na.rm = TRUE),
outlier = ifelse(meanIntensity > median + sd3 | meanIntensity < median - sd3, "yes", "no"),
dotProduct = ifelse(outlier == "yes", NA, dotProduct))
# plot after
dots.merged %>%
filter(algorithm == "logistic") %>%
ggplot(aes(1, meanIntensity)) +
geom_boxplot()trial.numbers = data.frame(subjectID = c(rep("DEV060", 79), rep("DEV061", 79), rep("DEV063", 71), rep("DEV081", 80), rep("DEV082", 75)),
rownum = c(1:79, 1:79, 1:71, 1:80, 1:75),
trial = c(1:19, 21:80, 1:59, 61:80, 1:31, 41:80, 1:20, 41:80, 21:40, 1:35, 41:80))
dots.check = dots.merged %>%
group_by(subjectID, algorithm) %>%
mutate(rownum = row_number()) %>%
left_join(., trial.numbers, by = c("subjectID", "rownum")) %>%
mutate(trial = ifelse(is.na(trial), rownum, trial),
dotProduct = ifelse(subjectID == "DEV060" & trial %in% 19:20, NA,
ifelse(subjectID == "DEV061" & trial == 59, NA,
ifelse(subjectID == "DEV082" & trial %in% 19:20, NA, dotProduct)))) %>%
select(-rownum) %>%
left_join(., striping, by = c("subjectID", "beta")) %>%
mutate(dotProduct = ifelse(!is.na(striping), NA, dotProduct))dots.ex = dots.check %>%
filter(!algorithm %in% c("ridge", "svm")) %>%
group_by(subjectID, algorithm) %>% # standardize within sub and algorithm
mutate(dotSTD = scale(dotProduct, center = FALSE)) Exclusions
Other * select only craved trials
data = left_join(dots.ex, data.ex, by = c("subjectID", "trial")) %>%
filter(!subjectID %in% c("DEV001","DEV020","DEV032","DEV047","DEV055","DEV064","DEV066", "DEV017", "DEV019", "DEV037", "DEV054")) %>%
filter(!(subjectID == "DEV029" & run == "run3") & !(subjectID == "DEV037" & run == "run1") & !(subjectID == "DEV042" & run == "run4") & !(subjectID == "DEV067" & run == "run4")) %>%
ungroup() %>%
mutate(algorithm = ifelse(algorithm == "reg_look", "regulate > look",
ifelse(algorithm == "reg", "regulate > rest", algorithm))) %>%
filter(craving == "craved")data %>%
filter(algorithm == "logistic") %>%
group_by(subjectID) %>%
summarize(n = n()) %>%
arrange(n)# roc curve
perf.df = data %>%
mutate(instruction = ifelse(instruction == "regulate", 1, 0)) %>%
group_by(algorithm) %>%
do({
algorithm = .$algorithm
pred = prediction(.$dotSTD, .$instruction, label.ordering = NULL)
perf = performance(pred, measure = "tpr", x.measure = "fpr")
data.frame(cut=perf@alpha.values[[1]],fpr=perf@x.values[[1]],tpr=perf@y.values[[1]])
})
ggplot(perf.df, aes(fpr, tpr, color = algorithm)) +
geom_line() +
scale_color_manual(values = algorithm) +
xlab("false positive rate") +
ylab("true positive rate")# logistic
roc = data %>%
filter(algorithm == "logistic") %>%
select(subjectID, trial, instruction, rating, dotSTD) %>%
mutate(guess.instruction = ifelse(dotSTD > 0, "regulate", "look"),
guess.rating = ifelse(dotSTD > 0, "low", "high"),
rating.bin = ifelse(rating >= 3, "high", "low"),
instruction = as.factor(instruction),
guess.instruction = as.factor(guess.instruction),
rating.bin = as.factor(rating.bin),
guess.rating = as.factor(guess.rating))
confusionMatrix(roc$guess.instruction, roc$instruction)## Confusion Matrix and Statistics
##
## Reference
## Prediction look regulate
## look 881 338
## regulate 529 1073
##
## Accuracy : 0.6927
## 95% CI : (0.6753, 0.7097)
## No Information Rate : 0.5002
## P-Value [Acc > NIR] : < 0.00000000000000022
##
## Kappa : 0.3853
##
## Mcnemar's Test P-Value : 0.0000000001098
##
## Sensitivity : 0.6248
## Specificity : 0.7605
## Pos Pred Value : 0.7227
## Neg Pred Value : 0.6698
## Prevalence : 0.4998
## Detection Rate : 0.3123
## Detection Prevalence : 0.4321
## Balanced Accuracy : 0.6926
##
## 'Positive' Class : look
##
confusionMatrix(roc$guess.rating, roc$rating.bin)## Confusion Matrix and Statistics
##
## Reference
## Prediction high low
## high 836 343
## low 768 750
##
## Accuracy : 0.5881
## 95% CI : (0.5692, 0.6067)
## No Information Rate : 0.5947
## P-Value [Acc > NIR] : 0.7661
##
## Kappa : 0.1953
##
## Mcnemar's Test P-Value : <0.0000000000000002
##
## Sensitivity : 0.5212
## Specificity : 0.6862
## Pos Pred Value : 0.7091
## Neg Pred Value : 0.4941
## Prevalence : 0.5947
## Detection Rate : 0.3100
## Detection Prevalence : 0.4372
## Balanced Accuracy : 0.6037
##
## 'Positive' Class : high
##
# regulate > look
roc = data %>%
filter(algorithm == "regulate > look") %>%
select(subjectID, trial, instruction, rating, dotSTD) %>%
mutate(guess.instruction = ifelse(dotSTD > 0, "regulate", "look"),
guess.rating = ifelse(dotSTD > 0, "low", "high"),
rating.bin = ifelse(rating >= 3, "high", "low"),
instruction = as.factor(instruction),
guess.instruction = as.factor(guess.instruction),
rating.bin = as.factor(rating.bin),
guess.rating = as.factor(guess.rating))
confusionMatrix(roc$guess.instruction, roc$instruction)## Confusion Matrix and Statistics
##
## Reference
## Prediction look regulate
## look 875 551
## regulate 535 860
##
## Accuracy : 0.615
## 95% CI : (0.5968, 0.633)
## No Information Rate : 0.5002
## P-Value [Acc > NIR] : <0.0000000000000002
##
## Kappa : 0.2301
##
## Mcnemar's Test P-Value : 0.649
##
## Sensitivity : 0.6206
## Specificity : 0.6095
## Pos Pred Value : 0.6136
## Neg Pred Value : 0.6165
## Prevalence : 0.4998
## Detection Rate : 0.3102
## Detection Prevalence : 0.5055
## Balanced Accuracy : 0.6150
##
## 'Positive' Class : look
##
confusionMatrix(roc$guess.rating, roc$rating.bin)## Confusion Matrix and Statistics
##
## Reference
## Prediction high low
## high 862 513
## low 742 580
##
## Accuracy : 0.5347
## 95% CI : (0.5156, 0.5536)
## No Information Rate : 0.5947
## P-Value [Acc > NIR] : 1
##
## Kappa : 0.0659
##
## Mcnemar's Test P-Value : 0.0000000001227
##
## Sensitivity : 0.5374
## Specificity : 0.5306
## Pos Pred Value : 0.6269
## Neg Pred Value : 0.4387
## Prevalence : 0.5947
## Detection Rate : 0.3196
## Detection Prevalence : 0.5098
## Balanced Accuracy : 0.5340
##
## 'Positive' Class : high
##
data %>%
filter(!is.na(rating)) %>%
ggplot(aes(algorithm, dotSTD, fill = instruction)) +
stat_summary(fun.y = mean, geom = "bar", position = position_dodge(width = 0.95)) +
stat_summary(fun.data = mean_cl_boot, geom = "errorbar", position = position_dodge(width = 0.95), width = 0) +
scale_fill_manual(name = "", values = instruction) +
labs(y = "standardized regulation pattern expression value\n", x = "") +
dc_bwdata %>%
filter(!is.na(rating)) %>%
ggplot(aes(instruction, dotSTD)) +
stat_summary(aes(group = subjectID), fun.y = mean, geom = "line", alpha = .1, size = .5) +
stat_summary(aes(group = 1), fun.y = mean, geom = "line", size = .75) +
stat_summary(aes(color = instruction), fun.data = "mean_cl_boot", geom = "pointrange", width = 0, size = .75) +
facet_grid(~algorithm) +
scale_color_manual(name = "", values = instruction) +
labs(y = "standardized regulation pattern expression value\n", x = "") +
dc_bwdata %>%
filter(!is.na(rating)) %>%
ggplot(aes(algorithm, dotSTD, fill = as.factor(rating))) +
stat_summary(fun.y = mean, geom = "bar", position = position_dodge(width = 0.95)) +
stat_summary(fun.data = mean_cl_boot, geom = "errorbar", position = position_dodge(width = 0.95), width = 0) +
scale_fill_manual(name = "", values = rating) +
facet_grid(~instruction) +
labs(y = "standardized regulation pattern expression value\n", x = "") +
dc_bwdata %>%
filter(!is.na(rating)) %>%
group_by(rating, algorithm, instruction) %>%
mutate(n.obs = n()) %>%
ggplot(aes(as.factor(rating), dotSTD, color = algorithm)) +
stat_summary(aes(group = algorithm), fun.y = mean, geom = "line") +
stat_summary(fun.data = "mean_cl_boot", geom = "linerange") +
stat_summary(aes(size = n.obs), fun.y = mean, geom = "point") +
scale_color_manual(name = "", values = algorithm) +
scale_size(name = "", range = c(1,4)) +
facet_grid(~instruction) +
labs(x = "\nrating", y = "standardized regulation pattern expression value\n") +
dc_bwdata %>%
filter(!is.na(rating)) %>%
group_by(rating, algorithm, instruction) %>%
mutate(n.obs = n()) %>%
ggplot(aes(as.factor(rating), dotSTD, color = instruction)) +
stat_summary(aes(group = instruction), fun.y = mean, geom = "line") +
stat_summary(fun.data = "mean_cl_boot", geom = "linerange") +
stat_summary(aes(size = n.obs), fun.y = mean, geom = "point") +
scale_color_manual(name = "", values = instruction) +
scale_size(name = "", range = c(1,4)) +
facet_grid(~algorithm) +
labs(x = "\nrating", y = "standardized regulation pattern expression value\n") +
dc_bwdata %>%
filter(!is.na(rating)) %>%
mutate(rating.bin = ifelse(rating >= 3, "high", "low")) %>%
ggplot(aes(rating.bin, dotSTD)) +
stat_summary(aes(group = subjectID), fun.y = mean, geom = "line", alpha = .1, size = .5) +
stat_summary(aes(group = 1), fun.y = mean, geom = "line", size = .75) +
stat_summary(aes(color = rating.bin), fun.data = mean_cl_boot, geom = "pointrange", width = 0, size = .75) +
facet_grid(~algorithm) +
scale_color_manual(name = "", values = craving) +
coord_cartesian(ylim = c(-1, 1.25)) +
labs(y = "standardized regulation pattern expression value\n", x = "") +
dc_bwdata %>%
filter(!is.na(rating)) %>%
mutate(rating.bin = ifelse(rating >= 3, "high", "low")) %>%
ggplot(aes(rating.bin, dotSTD)) +
stat_summary(aes(group = interaction(subjectID, instruction), color = instruction), fun.y = mean, geom = "line", alpha = .1, size = .5) +
stat_summary(aes(group = instruction, color = instruction), fun.y = mean, geom = "line", size = .75) +
stat_summary(aes(color = instruction), fun.data = mean_cl_boot, geom = "pointrange", width = 0, size = .75) +
facet_grid(~algorithm) +
scale_color_manual(name = "", values = instruction) +
coord_cartesian(ylim = c(-1.5, 2)) +
labs(y = "standardized regulation pattern expression value\n", x = "") +
dc_bwdata %>%
filter(!is.na(rating)) %>%
ggplot(aes(rt, dotSTD, color = instruction)) +
geom_smooth(method = "lm", alpha = .2) +
facet_grid(~algorithm) +
scale_color_manual(name = "", values = instruction) +
labs(x = "\nreaction time (seconds)", y = "standardized regulation pattern expression value\n") +
dc_bwdata %>%
filter(!is.na(rating)) %>%
filter(algorithm == "logistic") %>%
ggplot(aes(as.factor(rating), rt, color = instruction)) +
stat_summary(aes(group = instruction), fun.y = mean, geom = "line") +
stat_summary(fun.data = mean_cl_boot, geom = "pointrange", width = 0) +
facet_grid(~craving) +
scale_color_manual(values = instruction) +
labs(x = "\nrating", y = "reaction time (seconds)\n") +
dc_bwdata %>%
filter(!is.na(rating)) %>%
filter(algorithm == "logistic") %>%
ggplot(aes(rt, dotSTD, color = as.factor(rating))) +
geom_smooth(method = "lm", alpha = .1) +
facet_grid(~instruction) +
scale_color_manual(name = "rating", values = rating) +
labs(x = "\nreaction time (seconds)", y = "standardized regulation pattern expression value\n") +
dc_bwdata %>%
filter(!is.na(rating)) %>%
group_by(subjectID, algorithm, instruction, craving) %>%
mutate(meanPEV = mean(dotProduct, na.rm = TRUE)) %>%
ggplot(aes(instruction, meanPEV, color = instruction)) +
geom_boxplot() +
facet_wrap(~algorithm, scales = "free", ncol = 4) +
scale_color_manual(values = instruction) +
labs(x = "", y = "mean regulation pattern expression value\n") +
dc_bwdata %>%
filter(!is.na(rating)) %>%
group_by(subjectID, algorithm, instruction) %>%
mutate(meanPEV = mean(dotProduct, na.rm = TRUE),
meanRating = mean(rating, na.rm = TRUE)) %>%
select(subjectID, algorithm, meanPEV, meanRating, instruction) %>%
unique() %>%
ggplot(aes(meanPEV, meanRating, color = instruction)) +
geom_point() +
geom_smooth(method = "lm", se = .1) +
facet_wrap(~algorithm, scales = "free", ncol = 4) +
scale_color_manual(values = instruction) +
labs(x = "\nmean regulation pattern expression value", y = "mean craving rating\n") +
dc_bwdata %>%
filter(!is.na(rating)) %>%
group_by(subjectID, algorithm, instruction) %>%
mutate(meanPEV = mean(dotProduct, na.rm = TRUE),
meanRating = mean(rating, na.rm = TRUE),
instruction1 = instruction) %>%
select(subjectID, algorithm, meanPEV, meanRating, instruction, instruction1) %>%
unique() %>%
spread(instruction, meanRating) %>%
group_by(subjectID, algorithm) %>%
fill(everything(), .direction = "down") %>%
fill(everything(), .direction = "up") %>%
mutate(success = look - regulate,
success.percent = ((look - regulate) / look) * 100) %>%
spread(instruction1, meanPEV) %>%
mutate(diff = regulate - look) %>%
ggplot(aes(diff, success.percent)) +
geom_point() +
geom_smooth(method = "lm", se = .1) +
facet_wrap(~algorithm, scales = "free", ncol = 4) +
#scale_color_manual(values = instruction) +
labs(x = "\ndifference in mean regulation pattern expression value (regulate - look)", y = "percent change in craving (look - regulate / look)") +
dc_bwdata %>%
filter(!is.na(rating)) %>%
group_by(subjectID, algorithm, instruction) %>%
mutate(meanPEV = mean(dotProduct, na.rm = TRUE),
meanRating = mean(rating, na.rm = TRUE),
instruction1 = instruction) %>%
select(subjectID, algorithm, meanPEV, meanRating, instruction, instruction1) %>%
unique() %>%
spread(instruction, meanRating) %>%
group_by(subjectID, algorithm) %>%
fill(everything(), .direction = "down") %>%
fill(everything(), .direction = "up") %>%
mutate(success = look - regulate,
success.percent = ((look - regulate) / look) * 100) %>%
spread(instruction1, meanPEV) %>%
mutate(diff = regulate - look) %>%
ungroup() %>%
nest(-algorithm) %>%
mutate(
test = map(data, ~ cor.test(.x$success.percent, .x$diff)),
tidied = map(test, broom::tidy)
) %>%
unnest(tidied, .drop = TRUE)# load qualtrics
cred_file_location = '~/credentials.yaml.DEFAULT'
sid_column_name = 'Login|qid'
survey_name_filter = 'DEV Session 1 Surveys'
sid_pattern = 'DEV[0-9]{3}'
exclude_sid = '99|DEV737' # subject IDs to exclude
# load credential file
credentials = scorequaltrics::creds_from_file(cred_file_location)
# filter
surveysAvail = scorequaltrics::get_surveys(credentials)
surveysFiltered = filter(surveysAvail, grepl(survey_name_filter, SurveyName))
# get data
surveys_long = scorequaltrics::get_survey_data(surveysFiltered,
credentials,
pid_col = sid_column_name) %>%
mutate(Login = gsub("Dev", "DEV", Login),
Login = gsub("dev", "DEV", Login),
Login = ifelse(grepl("^[0-9]{3}$", Login), paste0("DEV", Login), Login),
Login = ifelse(Login == "DEVO55", "DEV055", Login))
# check for repeats
repeat.subs = surveys_long %>%
filter(item == "Finished" & value == "1") %>%
filter(!grepl(exclude_sid, Login) & !is.na(Login)) %>%
group_by(survey_name, Login) %>%
summarize(n = n()) %>%
filter(n > 1) %>%
spread(survey_name, n)
# surveys_long %>%
# filter(item == "StartDate") %>%
# filter(Login %in% repeat.subs$Login) %>%
# group_by(survey_name, Login) %>%
# mutate(n = n()) %>%
# filter(n > 1)
# filter out repeats
surveys = surveys_long %>%
filter(Login %in% unique(data$subjectID)) %>%
filter(!qid %in% c("R_1M01CRpEgQ9Sjzx", "R_3JfnLZ2XhekmvvF", "R_RUXzgKp7Sne7865")) %>%
select(-qid) %>%
rename("subjectID" = Login)
# get and score restraint scale
restraint = surveys %>%
filter(grepl("RS", item)) %>%
mutate(value = ifelse(value == "", NA, value),
value = as.integer(value)) %>%
group_by(subjectID) %>%
summarize(restraint = sum(value, na.rm = TRUE))
# correlation with other measures
data %>%
filter(!is.na(rating)) %>%
group_by(subjectID, algorithm, instruction) %>%
mutate(meanPEV = mean(dotProduct, na.rm = TRUE)) %>%
select(subjectID, algorithm, meanPEV, instruction) %>%
unique() %>%
group_by(subjectID) %>%
spread(instruction, meanPEV) %>%
mutate(diff = regulate - look) %>%
left_join(., restraint) %>%
ggplot(aes(diff, restraint)) +
geom_point() +
geom_smooth(method = "lm", se = .1) +
facet_wrap(~algorithm, scales = "free", ncol = 3) +
labs(x = "\ndifference in mean regulation pattern expression value (regulate - look)", y = "restraint score\n") +
dc_bwdata %>%
filter(!is.na(rating)) %>%
group_by(subjectID, algorithm, instruction) %>%
mutate(meanPEV = mean(dotProduct, na.rm = TRUE)) %>%
select(subjectID, algorithm, meanPEV, instruction) %>%
unique() %>%
group_by(subjectID) %>%
spread(instruction, meanPEV) %>%
mutate(diff = regulate - look) %>%
left_join(., restraint) %>%
ungroup() %>%
nest(-algorithm) %>%
mutate(
test = map(data, ~ cor.test(.x$restraint, .x$diff)),
tidied = map(test, broom::tidy)
) %>%
unnest(tidied, .drop = TRUE)cors = data %>%
filter(!is.na(rating)) %>%
group_by(subjectID, algorithm, instruction) %>%
mutate(meanPEV = mean(dotProduct, na.rm = TRUE),
meanRating = mean(rating, na.rm = TRUE)) %>%
gather(variable, value, meanPEV, meanRating) %>%
select(subjectID, instruction, algorithm, variable, value) %>%
unique() %>%
unite("instruction", c("variable", "instruction")) %>%
ungroup() %>%
select(subjectID, algorithm, instruction, value) %>%
group_by(algorithm) %>%
do({
instruction.spread = spread(., instruction, value)
cors = cor(instruction.spread[,-c(1:2)], use = "pairwise.complete.obs") %>%
as.data.frame() %>%
mutate(algorithm = instruction.spread$algorithm[[1]],
instruction = colnames(instruction.spread)[-c(1:2)])
})
cors %>%
reshape2::melt() %>%
ggplot(aes(instruction, variable, fill = value)) +
geom_tile(color = "white") +
scale_fill_gradientn(name = "correlation\n", colors = c("#3B9AB2", "white", "#F21A00"), limits = c(-1, 1), breaks = c(-1, 0, 1)) +
geom_text(aes(label = round(value, 2)), size = 3) + #family = "Futura Medium"
#geom_text(aes(label = significant), size = 4, family = "Futura Medium", nudge_x = .3, nudge_y = .2) +
facet_wrap(~algorithm) +
labs(x = "", y = "") +
dc_bw +
theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1))# set na.action for dredge
options(na.action = "na.fail")
# tidy data
data.sca = data %>%
select(subjectID, trial, algorithm, dotSTD, rating, instruction) %>%
spread(algorithm, dotSTD) %>%
na.omit()
# run full model
models = lmer(rating ~ instruction*logistic + instruction*`regulate > rest` + instruction*`regulate > look` + (1 | subjectID), data = data.sca)
# run all possible nested models
models.sca = dredge(models, rank = "AIC", extra = "BIC") %>%
select(AIC, delta, BIC, df, logLik, weight, `(Intercept)`, instruction, everything())
# set AIC for null model you want to compare model AIC values to
null = models.sca %>%
arrange(df) %>%
filter(instruction == "+" & df == 4)# tidy for plotting
plot.data = models.sca %>%
arrange(AIC) %>%
mutate(specification = row_number(),
better.fit = ifelse(AIC == null$AIC, "equal",
ifelse(AIC < null$AIC, "yes", "no")))
order = plot.data %>%
arrange(AIC) %>%
mutate(better.fit.num = ifelse(better.fit == "yes", 1, 0)) %>%
gather(variable, value, -c(AIC, delta, BIC, df, logLik, weight, better.fit.num, specification, better.fit)) %>%
filter(!is.na(value)) %>%
group_by(variable) %>%
mutate(order = sum(better.fit.num)) %>%
select(variable, order) %>%
unique()
# variables included in model
variable.names = names(select(plot.data, -starts_with("better"), -specification, -AIC, -BIC, -df, -logLik, -delta, -weight, -`(Intercept)`))
# plot top panel
a = plot.data %>%
ggplot(aes(specification, AIC, color = better.fit)) +
geom_point(shape = "|", size = 4, alpha = .75) +
geom_hline(yintercept = null$AIC, linetype = "dashed", color = "#5BBCD6") +
scale_color_manual(values = c("#5BBCD6", "black", "#F43C13")) +
labs(x = "", y = "AIC\n") +
theme_minimal(base_size = 14) +
theme(legend.title = element_text(size = 10),
legend.text = element_text(size = 9),
axis.text = element_text(color = "black"),
axis.line = element_line(colour = "black"),
legend.position = "none",
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.border = element_blank(),
panel.background = element_blank())
# set plotting order for variables based on number of times it's included in better fitting models
b = plot.data %>%
gather(variable, value, eval(variable.names)) %>%
left_join(., order, by = "variable") %>%
mutate(value = ifelse(!is.na(value), "|", ""),
variable = gsub("`regulate > look`", "regulate > look", variable),
variable = gsub("`regulate > rest`", "regulate > rest", variable),
variable = gsub("regulate > look:instruction", "instruction:regulate > look", variable),
variable = gsub("regulate > rest:instruction", "instruction:regulate > rest", variable),
variable = gsub("instruction:", "instruction x ", variable)) %>%
ggplot(aes(specification, reorder(variable, order), color = better.fit)) +
geom_text(aes(label = value), alpha = .75) +
scale_color_manual(values = c("#5BBCD6", "black", "#F43C13")) +
labs(x = "\nspecification number", y = "variables\n") +
theme_minimal(base_size = 14) +
theme(legend.title = element_text(size = 10),
legend.text = element_text(size = 9),
axis.text = element_text(color = "black"),
axis.line = element_line(colour = "black"),
legend.position = "none",
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.border = element_blank(),
panel.background = element_blank())
(plot = plot_grid(a, b, ncol = 1, align = "v"))